! ---
! Copyright (C) 1996-2016	The SIESTA group
!  This file is distributed under the terms of the
!  GNU General Public License: see COPYING in the top directory
!  or http://www.gnu.org/copyleft/gpl.txt .
! See Docs/Contributors.txt for a list of contributors.
! ---
      MODULE m_state_analysis
      use write_subs

      private
      public :: state_analysis

      CONTAINS

      subroutine state_analysis( istep )
      use siesta_cml
      use m_born_charge, only : born_charge
      use parallel,      only : IOnode
      use m_wallclock,   only : wallclock
      use zmatrix,       only : lUseZmatrix, iofaZmat,
     &                          CartesianForce_to_ZmatForce
      use atomlist,      only : iaorb, iphorb, amass, no_u, lasto
      use atomlist,      only : indxuo, Qtot
      use m_spin,        only : spin
      use m_fixed,       only : fixed
      use sparse_matrices
      use siesta_geom

      USE siesta_options
      use units, only: amu, eV
      use m_stress
      use m_energies, only: Etot, FreeE, Eharrs, FreeEHarris, Entropy
      use m_energies, only: Ebs, Ef
      use m_ntm
      use m_forces
      use m_energies, only: update_FreeE, update_FreeEHarris
      use m_intramol_pressure, only: remove_intramol_pressure
      use m_ts_global_vars, only: TSrun
      use m_ts_options, only: N_elec, Elecs
      use ts_charge_m, only: ts_charge_print, TS_Q_INFO_FULL

#ifdef SIESTA__FLOOK
      use flook_siesta, only : slua_call, LUA_FORCES
#endif

      implicit none

      integer  :: istep
      integer  :: ia, jx, ix
      real(dp) :: volume
      logical  :: eggbox_block=.true. ! Read eggbox info from data file?
      real(dp) :: qspin
      
      external :: eggbox, mulliken, moments
      real(dp), external :: volcel

!------------------------------------------------------------------------- BEGIN
      call timer( 'state_analysis', 1 )
#ifdef DEBUG
      call write_debug( '  PRE state_analysis' )
#endif

      if (cml_p) then
        call cmlStartModule(xf=mainXML, title='SCF Finalization')   
      endif

!     Write final Kohn-Sham and Free Energy
      
      FreeE       = Etot - Temp * Entropy
      FreeEHarris = Eharrs - Temp * Entropy

      if (cml_p) call cmlStartPropertyList(mainXML,
     &                        title='Energies and spin')
      if (IOnode) then
        if ( .not. harrisfun)
     &      write(6,"(/a,f14.4)")  'siesta: E_KS(eV) =        ', Etot/eV
        if (cml_p) then
           call cmlAddProperty(xf=mainXML, value=Etot/eV,
     &       dictref='siesta:E_KS', units='siestaUnits:eV', 
     .       fmt='r6')
           call cmlAddProperty(xf=mainXML, value=FreeE/eV,
     &       dictref='siesta:FreeE', units='siestaUnits:eV', 
     .       fmt='r6')
           call cmlAddProperty(xf=mainXML, value=Ebs/eV,
     &       dictref='siesta:Ebs', units='siestaUnits:eV', 
     .       fmt='r6')
           call cmlAddProperty(xf=mainXML, value=Ef/eV,
     &       dictref='siesta:E_Fermi', units='siestaUnits:eV', 
     .       fmt='r6')
        endif
      endif

!     Substract egg box effect from energy
      if (eggbox_block) then
        call eggbox( 'energy', ucell, na_u, isa, ntm, xa, fa, Etot,
     &               eggbox_block )
        FreeE  = Etot - Temp * Entropy
        if (IOnode)
     &    write(6,"(/a,f14.4)") 'siesta: E_KS - E_eggbox = ',Etot/eV
        if (cml_p) call cmlAddProperty(xf=mainXML, value=Etot/eV,
     &         dictref='siesta:E_KS_egg', units='siestaUnits:eV', 
     .         fmt='r6')
      endif

      call update_FreeE( Temp )
      call update_FreeEHarris( Temp )
      call print_spin(qspin)
      
      if (cml_p) call cmlEndPropertyList( mainXML )

!     Substract egg box effect from the forces 
      if (eggbox_block) then
        call eggbox('forces',ucell,na_u,isa,ntm,xa,fa,Etot,eggbox_block)
      endif

      if (IOnode) call write_raw_efs(stress,na_u,fa,FreeE)

!     Compute stress without internal molecular pressure
      call remove_intramol_pressure(ucell,stress,na_u,xa,fa,mstress)

!     Impose constraints to atomic movements by changing forces ...........
      if (RemoveIntraMolecularPressure) then
!        Consider intramolecular pressure-removal as another
!        kind of constraint
         call fixed( ucell, mstress, na_u, isa, amass, xa, fa,
     &               cstress, cfa, ntcon , 
     &               magnitude_usage = idyn==0 )
      else
         call fixed( ucell, stress, na_u, isa, amass, xa, fa,
     &               cstress, cfa, ntcon ,
     &               magnitude_usage = idyn==0 )
      endif

#ifdef SIESTA__FLOOK
      ! We call it right after using the
      ! geometry constraints.
      ! In that way we can use both methods on top
      ! of each other!
      ! The easy, already implemented methods in fixed,
      ! and custom ones in Lua :)
      call slua_call(LUA, LUA_FORCES)
#endif

!     Calculate and output Zmatrix forces
      if (lUseZmatrix .and. (idyn.eq.0)) then
        call CartesianForce_to_ZmatForce(na_u,xa,fa)
        if (IOnode) call iofaZmat()
      endif

!     Compute kinetic contribution to stress
      kin_stress(1:3,1:3) = 0.0_dp
      volume = volcel(ucell)
      do ia = 1,na_u
        do jx = 1,3
          do ix = 1,3
            kin_stress(ix,jx) = kin_stress(ix,jx) -
     &             amu * amass(ia) * va(ix,ia) * va(jx,ia) / volume
          enddo
        enddo
      enddo
!     Add kinetic term to stress tensor
      tstress = stress + kin_stress

!     Force output 
      if (IOnode) then
        call siesta_write_forces(istep)
        if ( TSrun ) then
          call transiesta_write_forces()
        end if
        call siesta_write_stress_pressure()
        call wallclock('--- end of geometry step')
      endif

!     Population and moment analysis 
      if ( spin%SO .and. orbmoms) then
         call moments( 1, na_u, no_u, maxnh, numh, listhptr,
     .           listh, S, Dscf, isa, lasto, iaorb, iphorb,
     .           indxuo )
      endif
      ! Call this unconditionally
      call mulliken( mullipop, na_u, no_u, maxnh,
     &               numh, listhptr, listh, S, Dscf, isa,
     &               lasto, iaorb, iphorb )

      ! Also write TS charges
      if ( TSrun ) then
        call ts_charge_print(N_Elec,Elecs, Qtot,
     &      block_dist, sparse_pattern,
     &      spin%DM, maxnh, Dscf, S, TS_Q_INFO_FULL)
      end if

!
!     Call the born effective charge routine only in those steps (even) 
!     in which the dx  is positive.
      if (bornz .and. (mod(istep,2) .eq. 0)) then
         call born_charge()
      endif

!     End the xml module corresponding to the analysis
      if (cml_p) then
         call cmlEndModule(mainXML)         
      endif 
      call timer( 'state_analysis', 2 )

!--------------------------------------------------------------------------- END
      END subroutine state_analysis

      END MODULE m_state_analysis
